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ABSTRACT 

Context. An Equation of State (EoS) is a relation between thermodynamic state variables and it is essential for closing the set 
of equations describing a fluid system. Although an ideal EoS with a constant adiabatic index F is the preferred choice due to its 
simplistic implementation, many astrophysical fluid simulations may benefit from a more sophisticated treatment that can account for 
diverse chemical processes. 

Aims. In the present work we first review the basic thermodynamic principles of a gas mixture in terms of its thermal and caloric EoS 
by including effects like ionization, dissociation as well as temperature dependent degrees of freedom such as molecular vibrations 
and rotations. The formulation is revisited in the context of plasmas that are either in equilibrium conditions (local thermodynamic- 
or collisional excitation- equilibria) or described by non-equilibrium chemistry coupled to optically thin radiative cooling. 

We then present a numerical implementation of thermally ideal gases obeying a more general caloric EoS with non-constant adiabatic 
index in Godunov-type numerical schemes. 

Methods. We discuss the necessary modifications to the Riemann solver and to the conversion between total energy and pressure (or 
vice-versa) routinely invoked in Godunov-type schemes. We then present two different approaches for computing the EoS. The first 
one employs root-finder methods and it is best suited for EoS in analytical form. The second one leans on lookup table and interpola¬ 
tion and results in a more computationally efficient approach although care must be taken to ensure thermodynamic consistency. 
Results. A number of selected benchmarks demonstrate that the employment of a non-ideal EoS can lead to important differences in 
the solution when the temperature range is 500 - 10“* K where dissociation and ionization occur. The implementation of selected EoS 
introduces additional computational costs although the employment of lookup table methods (when possible) can significantly reduce 
the overhead by a factor 3 ~ 4. 
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1. Introduction 


An equation of state (EoS) is a relationship between state vari¬ 
ables of a thermodynamic system under certain physical condi¬ 
tions. Such a constitutive equation provide a necessary closure 
for a complete mathematical description of a fluid system in 
addition to the conservation laws of mass, momentum and en¬ 
ergy. Numerical simulations of astrophysical systems such as 
inter-stellar medium, planetary atmospheres, stellar evolution, 
jets and outflows, require inter-play of various thermal, radia¬ 
tive and chemical processes. Eor such complex systems, using 
a simple ideal (or an isothermal) EoS would be considered as a 
serious limitation. A consistent description for such systems de¬ 
mands the use of a general EoS that can account for thermal and 
chemical processes. 


Eor example, the thermodynamic state of the gas plays a piv¬ 
otal role in governing the fragmentation of self-gravitating and 
turbulent molecular clouds (e.g., Spaans & Silk 2000 [ Li et al. 


2003 1 Jappsen et al. 2005|l. The balance between heating and 


cooling in molecular clouds is approximated by using a poly¬ 
tropic EoS, p oc . Multiple smoothed particle hydrodynami- 
cal simulations with different adiabatic indices, 0.2 < E < 1.4 
( |Spaans & Silk|2000| ) were used to show that the degree of frag¬ 
mentation decreases with increasing value of E ( |Li et al. 12003) 1. 
Jappsen et al. ( 2005|l showed that the thermal properties of the 


gas determines the stellar mass function (IMF) using a piece- 


wise poly-tropic EoS. Such empirical forms of EoS in general 
depend on chemical abundances and complex atomic and molec¬ 
ular physics. 

Numerical simulations studying thermo-chemical evolution 
of early structure formation used an effective adiabatic index. 


Eeff, to relate internal energy with thermal pressure (e.g. Yoshida 
|et al.|20()6l [Glover & Abel|2008| l. The value of Eeff is estimated 
from number fractions of chemical species treating the chemical 
composition as an ideal mixture. In the context of disk instabil¬ 
ity leading to formation of gas giant planets, [Boley et aT (|2007 1 


pointed out the importance of incorporating isotopic forms of 
molecular hydrogen, H 2 , as well the molecular physics (rota¬ 
tion and vibration) under thermodynamic equilibrium in the es¬ 
timate of internal energy. A more complex EoS taking into ac¬ 
count ionization of atomic hydrogen, helium and radiation along 
with molecular dissociations is used to study the envelopes of 
young planetary cores (D’Angelo & Bodenheimer|2013|l. 


From the computational perspective, procedures to treat dy¬ 
namics of astrophysical plasmas with a general EoS have been 


described, e.g., by 

Colella & GlazI (fT985ll; Glaister (1988b|a|l; 

Menikoff & Plohr 

(|1989|l; |Liou et al.| (|1990|l; |Fedkiw et al.| 

(1997i; 

Hu et al. (20091. Similarly, numerical codes like FLASH 

(Fryxel 

et al.pOOO)! and CASTRO (Almgren et al.||2010 

1 have 


implemented electron-positron EoS (|Timmes & Swesty 


20001 


using high order polynomials as interpolating functions. Such 
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an EoS based on tabulated Helmholtz free energy is employed 
in the study of stellar evolution and supernovae. Falle & Raga 
( 1993| 1995| l incorporated LTE effects to model partially ion¬ 
ized hydrogen gas using Saha Equations to study variable knots 
produced in stellar jets. 

The goal of this paper is to outline a consistent numerical 
framework for the implementation of a more general equation of 
state in the context of the magnetohydrodynamics (MHD) equa¬ 
tions. Our formulation accounts for different physical processes 
such as atomic ionization and recombination, molecular dissoci¬ 
ation, etc... and it is suitable under equilibrium conditions (local 
thermodynamic or collisional ionization equilibria) as well as 
for non-equilibrium optically thin radiative cooling ( |Te§ileanu 
et al.|200^ . The numerical method is implemented as part of the 
PLUTO code ( Mignone et al.|[2007| l and it is built while ensur¬ 
ing thermodynamical consistency, accuracy and computational 
efficiency. 

Our starting point are the ideal MHD equations written in 
conservation form: 


dp 

— -H V ■ 
dt 


(pv) = 0 


( 1 ) 


djpv) 

dt 


+ V ■ (pvv^ - + Vp, = 0 


( 2 ) 


^-Vx(vxB) = 0 

dt 


(3) 


a relation among intensive and extensive variables. A thermal 
EoS is an expression relating pressure, temperature and volume 
and we will express it as p = p{V, T). Conversely, a caloric EoS 
specifies the dependence of the internal energy of the system U 
on volume V and temperature T. The total internal energy, U is 
related to the internal energy density (see Eq.|^ as UjV - pe. 

In general, different forms of EoS relations are derived from 
empirical results and are used to estimate various thermodynam¬ 
ical properties of a system. Theoretically, statistical principles 
can be applied to describe such a system on basis of its micro¬ 
scopic processes using the partition function X- For example, the 
macroscopic thermodynamic quantities can be obtained from the 
following standard relations. 


p = ksT 


ld\nZ\ 

I j. 


U = kuT^ 


<91nZ 

dT 


( 8 ) 


where, ks is the Boltzmann constant. The previous equation es¬ 
sentially provides two forms of EoS in terms of partition func¬ 
tion. 

An important quantity that relates the pressure p with density 
p is the speed of sound, defined as 


Cs = 



(9) 


dE 

—+V-[(£ + p,)v-(vB)B] = A 
dt 


dipX,) 

dt 


+ V ■ ipXjpi) - S A 


(4) 


(5) 


where p is the mass density, v is the fluid velocity, B is the mag¬ 
netic field, Pt - P + B^/2 is the total pressure accounting for 
thermal (p) and magnetic (B^/2) contributions. The total energy 
density E is given by 


1 , 1 , 

E^pe+ -p\^ + -B^ . 


(6) 


An additional EoS relating the internal energy density pe with 
p and p must be specified. This issue is addressed in ^ 
Dissipative effects have been neglected for the sake of exposi¬ 
tion although they can be easily incorporated in this framework. 

The paper is organized as follows, in ^the basic principles 
and formulations of general EoS used for the present work are 
described. The numerical framework is discussed in ^ The re¬ 
sults obtained from various test problems are outlined in ^and 
the concluding remarks are summarized in ^ . 


2. Equation of State 

2.1. Thermodynamical Principles 

The principle of conservation of energy in thermodynamics is 
commonly known as the first law of thermodynamics and can be 
expressed as, 

dU^TdS-pdV, (7) 

where S(U,V) is the entropy. The internal energy U and the vol¬ 
ume V are classified as extensive variables and depend on bulk 
properties of the system. Whereas, the intensive variables like 
temperature T and pressure p show no dependence on the size 
of the system. An EoS describing such a system is defined as 


where the derivative must be taken at constant entropy. The 
above definition can be further expressed in terms of the first 
adiabatic exponent, Ei, defined as (51np/cHnp)s. Thus, Eq. 
now becomes 


( 10 ) 


In general the first adiabatic exponent Ei has a functional depen¬ 
dence on temperature and density as, 

1 ip\ 2^ 

■' \Xt 



r, = 


( 11 ) 


Cv(T)\pTi 

where, Cv{T) is obtained by taking the derivative of the specific 
gas internal energy, e(T) with respect to temperature at constant 
volume while andxp are referred to as temperature and den¬ 
sity exponents (see [D’Angelo & Bodenheimer ( |2013 i) 

/51np\ ^ cHnp 

“ l^lnT/^ “ " 51nT 

( 12 ) 

/cHnp\ ^ cHnp 
\(91npj^ dlnp' 

For an ideal gas, the value of Ei coincides with the adiabatic 
index E, which is essentially the ratio of specific heats. In such a 
case, the sound speed can be expressed as 


Xt 


Xp 


c, = 



(13) 


In the present work, we will focus on thermally ideal gases. 
These gases have their thermal EoS same as that of an ideal 
gas. However, the caloric EoS can have nonlinear dependence 
on temperature based on various chemical processes taken into 
consideration (see §2.3 i. We point out that, although the analy¬ 
sis presented here is limited to thermally ideal gas, the numerical 
implementation described in this work can also be extended to 
describe real gases obeying EoS that are not thermally ideal. 
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2.2. Thermodynamic Constraints 

The equation of state must adhere to a number of physical prin¬ 
ciples in order to be thermodynamically consistent. Although a 
comprehensive discussion lies outside the scope of this paper, we 
briefly outline the most relevant ones ( |Menikoff & Plohr|1989| l. 

1. The specific internal energy as a function of specific volume 
and entropy e - e(V,S) must be piecewise twice continu¬ 
ously differentiable. 

2. Thermodynamic stability demands that e(V,S) be a jointly 
convex function. This implies that the Hessian matrix of sec¬ 
ond derivatives of e with respect to V and S is non-negative. 

3. Simple physical considerations lead to the following asymp¬ 
totic conditions; 

lim p(V,S) = 0 (14) 


limpCVi^) = lim £>(1/5)= lim e(y, s) = oo (15) 

V—>0 5^00 5-^00 

The previous constraints also guarantee the existence of the 
solution of the Riemann problem. 


An additional constraint, convexity, can be introduced if one 
sticks to standard theory and phase transition are not considered. 
The convexity of an EoS is quantified by the fundamental gas 
derivative which expresses the nonlinear variation of the sound 
speed with respect to density and it is denoted by Q-. 


l + -(^ 

Cs \ op 


(16) 


In Eq. ( fTfll l the derivative is taken at constant entropy and Cj is 
the speed of sound given by Eq. (j^. Eor an ideal polytropic gas, 
one finds ^ = (T -H l)/2 while a general expression in terms of 
derivatives with respect to temperature and density may be found 
in Appendix [A| 

An EoS is said to be convex if ^ > 0 and it has the follow¬ 
ing important implications; i) the isoentropes are convex func¬ 
tions in the p -V plane, ii) the sound speed increases with den¬ 
sity along isoentropes, iii) only regular waves (e.g. compression 
shock waves and expansion fans) can be formed in the Riemann 
problem. 

The latter property is of particular interest here since, as we 
shall see, inaccurate table interpolation can lead to local viola¬ 
tion of the convexity assumption (see the discussion in Section 
3.2.3 and the results in Appendix [B]l. In such cases, compos¬ 


ite (or compound) waves consisting of a rarefaction wave prop¬ 
agating adjacent to a shock may be generated in the solution 
whilst satisfying the thermodynamical principles (e.g. Menikoff 
|& Plohr||l989| and references therein). This circumstance may 
arise, for instance, for a real gas in correspondence of finite in¬ 
tervals of concave p - V isoentropes. 

In the present paper, however, we restrict our attention to 
EoS for which ^ > 0 is always verified at the continuous level al¬ 
though it may not be true at the discrete numerical level thereby 
generating spurious composite waves. We address this issue in 
Sect. 1223] 


2.3. Caioricaiiy ideai Gas 


Consider the case of a classical monoatomic ideal gas, where, 
the partition function iZ is given by 



2nh^ j 



(17) 


where, m is the mass of the particle, h the Planck constant and 
N the total number of non-interacting particles. On substituting 
Eq. (17 1 in Eq. ([^, we obtain the standard EoS for a classical 
ideal gas, 

p = nksT 


pe = UjV = ^^nkeT, 


(18) 


where, n — N/V is the number density and /tj^ns denotes the 
translational degree of freedom which for a monoatomic gas 
equals to 3. Eurther, the specific heat capacity at constant vol¬ 
ume, Cv - /tians^/2 (where R being the universal gas constant), 
is independent of the temperature. On extending this analysis 
further to diatomic ideal gas, the partition function Z contains 
contribution from rotational and vibrational degrees of freedom, 
in addition to the translational motion. In such a case, the inter¬ 
nal energy density derived from Eqj^is given by. 


pe = ^nksT + ^nkgT + (Dv.b(r) (19) 


where the additional contribution of foinksT 12 comes from 
/rot rotational degree of freedoms, whose value is 2 for linear 
molecules and 3 for non-linear ones. In addition, <I>vib(T) de¬ 
notes term due to vibrational motion which has a non-linear de¬ 
pendence on temperature. On considering the diatomic molecule 
with two degrees of freedom (i.e, translational and rotational) 
and neglecting the non-linear dependence due to vibration, one 
obtains a single relation for both monoatomic and diatomic gas 
by adopting a constant E, 


p = (r-l)pe, (20) 

where T = 5/3 for monoatomic gas and E = 7/5 for diatomic 
gas (see Eqns. [TS] and [T^ . 


2.4. Partiaiiy ionized Hydrogen Gas 

Astrophysical fluids and processes are more complex than the 
simple system of ideal gas described above. Eor example, the 
Inter-Stellar Medium (ISM) that largely comprises hydrogen and 
helium is affected by many physical and chemical processes 
viz., collisional ionization, dissociation, shocks, radiation, etc. 
In such a scenario, an heuristic approach that models the ISM 
as an monoatomic ideal gas with constant E = 5/3 will only be 
approximate and fail to account for the feedback of the above 
processes on the thermal properties of the gas and thereby also 
on its inter-linked dynamics. 

Consider the simplest case of a partially ionized gas of pure 
hydrogen (in atomic form). The thermodynamics of such a sys¬ 
tem is different from that of a completely ionized (or completely 
neutral) gas as the number of free particles can change and an 
additional energy contribution is required during the process 
of ionization. The internal energy density is therefore given by 
( |Clayton|1984l l, 

pe = ^nkgT + Xunmi ■ (21) 

In addition to the standard form of translational energy, contribu¬ 
tion from ionization potential,^//, is included in Eq. ( |2T] l. Here, 
Whh is the number density of ionized hydrogens and the total 
number density of free particles, n - n^ + 2n////, is the sum of 
number densities of neutral hydrogens and twice that of nnn due 
to charge neutrality. 
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In regions of dense stellar interior, one can assume local ther¬ 
modynamic equilibrium. For such a system, the fractions of ion¬ 
ized hydrogen becomes nonlinearly dependent on temperature 
and density of the gas through the Saha equation. As a result, 
internal energy, specihc heats and the adiabatic index F will de¬ 
pend on the ionization fraction. For example, the adiabatic in¬ 
dex F will smoothly change from its monoatomic value of 5/3 
to 1.13 typical of a hydrogen gas with an ionization fraction of 
50% at r = 10"^ K ( Clayton|l984| . Such a signihcant change in 
F owes to the fact that part of the energy input becomes avail¬ 
able to ionization rather than increasing the temperature of the 
gas. Therefore, using a constant value of F for such dense stellar 
interiors will considerably overestimate the temperature of the 
gas. 


2.5. Hydrogen/Helium gas mixture 


with a mass fraction of Y and negligible fraction of metals. 
For such a composition the total density of gas is dehned as 
p - npifiH, where the mean molecular weight p can be expressed 
as (e.g., [Black & Bodenheimer|1975| l 

'^^[2X(\+y + 2xy) + Y]-\ (26) 

Such a gas mixture is further assumed to be thermally ideal so 
that pressure and temperature are related hy p - pksTKpiriH). 

The most crucial part is to express a caloric EoS that can 
account for contributions from various degrees of freedom and 
processes like ionization and dissociation. Thus, the gas internal 
energy density (pe)gas for the mixture is given by 

pksT 

(pe)gas - (pm + Chi + shu + sh+h + shc) -, (27) 

niH 


In recent years, studies related to planet formation in accretion 
disks have started to incorporate EoS that can account for con¬ 
tributions from dissociation of molecular hydrogen, ionization 


of atomic hydrogen and helium and radiation (e.g., Boley et al. 
2007 1 D’Angelo & Bodenheimer|2013 i under an assumption of 


local thermodynamic equilibrium (LTE). In this pwer, we have 
implemented such an EoS both in presence of LTI^and with ex¬ 
plicit non-equilibrium cooling. For all future references to this 
EoS, we will use H/He EoS. 

In LTE, processes like ionization-recombination and 
dissociation-bond formation for hydrogen are given by. 


H + e^ + 2e^ 

H2^H + H, 


( 22 ) 


respectively. Eollowing D’Angelo & Bodenheimer ( 2013| l, we 
define the degree of dissociation y and degree of ionization x as. 


Phi 

Phi + Ph2 

PHII 

Phi + PHI! ’ 


(23) 


where, pni is the density of atomic hydrogen, the density of 
molecular hydrogen and pun the density of ionized hydrogen. 
In the limit of LTE, one assumes that the level populations due 
to ionization (and dissociation) processes follow Boltzmann ex¬ 
citation formula and that the ejected free electrons thermalize to 
attain a Maxwell-Boltzmann velocity distribution corresponding 
to single gas temperature. This is generally true in regions of 
high density like that of the solar interior. In such cases, the de¬ 
gree of ionization using Saha equations is given as follows. 


_ ///«_ / ff^ekBT \ ^ B.eOcV/C/tflD 

1 - X Xp \ 2Kh^ ) 


(24) 


and also degree of dissociation, y can be obtained in a similar 
manner ([Black & Bodenheimer|1975|), 


_ ^H I niuksE X ^ AAieVKheT) 

I — y 2Xp \ Anh? ) 


(25) 


The gas is essentially a mixture of hydrogen in all forms 
(atoms, ions & molecules) with a mass fraction of X, Helium 


' For the present work, we have not considered contributions from 
radiation and He ionization but rather focused on hydrogen alone to 
be consistent with the implementation in presence of cooling which as¬ 
sumes a fixed fraction of helium. 


where each term in parenthesis is dimensionless and can be ob¬ 
tained from an appropriate partition function Z, and Eq. (18i. 
Table summarizes the different contribution to the gas internal 
energy. 

In the case of molecular hydrogen, eu,, terms that correspond 
to vibrational and rotational degree of freedom are also consid¬ 
ered. These terms are evaluated using the partition function of 
vibration and rotation that have explicit and a non-linear de¬ 
pendence on temperature. Additionally, the rotational partition 
function also takes into account the para/ortho H 2 spin states 
( Boley et al.[|2007] l. Thus, the total gas internal energy density 
has a nonlinear dependence on the temperature T and density 
through X and y (see Eqns. 24 &[25|). 

The left and middle panels of Eig. show the variation of 
pip, T) (Eq. [ 26) 1 and gas internal energy in ergs with tempera¬ 
ture, T, for four values of density in gcm~^ respectively. The 
values of p are bounded between the upper value ~2.3, corre¬ 
sponding to a fully molecular medium at low temperatures and a 
lower value ~ 0.6 at high temperatures representing a fully ion¬ 
ized medium. The transition between these bounds is smooth at 
large densities p - 10“"^ gcm“^ while it forms an intermediate 
plateau at T ~ 10^ K at low density values (black curve). The 
first transition occurs in the temperature range where molecules 
begin to dissociate to form atomic hydrogen. A second transi¬ 
tion takes place where atomic hydrogen becomes ionized. The 
same transitions can be observed in the prohle of internal en¬ 
ergy. Erom a physical point of view they indicate that the energy 
at these temperatures becomes available to dissociate or ionize 
the gas rather than heating the gas so that temperature remains 
approximately constant. Away from these transition regions, the 
dependence of (pe)gas(OT/r/p) is linear and increases monotoni- 
cally with the gas temperature. The last panel of the same hgure 
shows the variation of first adiabatic exponent, F 1 with temper¬ 
ature. At low temperatures, the gas behaves as a monoatomic 
ideal gas undergoing adiabatic process with Fi = 5/3. This is 
also true at very large temperatures where the gas contains ions 
and electrons. From the previous considerations, we see a sharp 
decrease in Fi from its maximum value of 5/3 to values around 
unity (corresponding to an isothermal limit) for a low density 
plasma (black curve). On the other hand, a single dip at T > 10"^ 
K is seen at larger densities (red curve). 

In addition to the study of planet formation in accretion 
discs, the H-He EoS is an important ingredient in the physics of 
proto-stellar formation from collapse of dense molecular cores. 
Radiation hydrodynamics simulations ( [Masunaga et al. 1998t 
Masunaga & Inutsuka|2000 1 have put forth detailed understand¬ 


ing of thermodynamics in presence of gravitational collapse. At 
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Fig. 1. Variation of Mean molecular weight yu, internal energy density of the gas (pe)gas and first adiabatic index T\ with temperature. The different 
colored curves represent four values of fixed density in gcm“^, viz., 10“'* (red),10“* (green), 10“*^ (blue) and iQ^^^(black). The values of (pe)gas 
and Ti are obtained at equilibrium between ortho and para hydrogen. 


Table 1. Summary of different contributions to the gas internal energy (pe)gas 
( |1975^ ; [D’Angelo & Bodenheimer| ( |2013^ ) 


, which is expressed using Eq. 1 27 


(see 


Black & Bodenheimer 


Term 

Expression 

Description 

eui 

1.5X(1 -Hx)y 

Translational energy for hydrogen 

SHu 

0.375Y 

Translational energy for helium 

rn+H 

4.48 tYXyKlksT) 

Dissociation energy for molecular hydrogen 

£hh 

13.6eYXxyl(ksT) 

Ionization energy for atomic hydrogen 

£h2 

X(l-a) \] IZ . T dp T dCr) 

2 V Cv dT ^ Cr dT \ 

Internal energy for molecular hydrogen 


the onset of collapse, compressional heating in the dense (p ~ 
10"^“^ g cm“^) and cold (T ~ lOK) core increases the central 
temperature up to T ~ 100 K adiabatically. As the temperature 
increases further, rotational states of H 2 are excited and the sys¬ 
tem evolves with an effective adiabatic index of 7/5 with the en¬ 
suing formation of a pressure-supported first core. Further in¬ 
crease of temperature beyond 10^ K results in dissociation of H 2 
which acts as an efficient cooling mechanism leading to a second 
collapse. 

However, not all astrophysical problems can be treated in 
LTE limit. A classical case is that of a jet, where the recombi¬ 
nation time scales are comparable to that of dynamical time. In 
such a scenario, LTE assumptions become invalid and a non¬ 
equilibrium approach has to be adopted as described in the fol¬ 
lowing section. 

2.6. Non-Equilibrium Hydrogen Chemistry 

Astrophysical flows in HII regions, supernova remnants, star 
forming regions are some classical examples where optically 
thin cooling time scales are comparable to the dynamical time. 
In such environments, ionization and dissociation fractions are 
far from LTE and their estimation based on Saha fractions can 
give large errors. In such cases, the number density of various 
species is more accurately determined by solving the chemical 
rate equations: 


where n is the number density, 'Kj^k is the rate of formation of t'* 
specie from all j and k species while TC, y is the rate of destruction 
of the f* specie due to all j species. 

In addition, proper treatment should be carried out to evolve 
the internal energy to account for losses due to optically thin 
radiation: 


d(pe) 

dt 


-A(n,X, T), 


(29) 


where A(n, X, T) is the optically thin radiative loss term. 
Radiative losses imply that the emitted photons due to differ¬ 
ent physical processes (e.g., ionization, metal line cooling etc.) 
freely stream (without diffusion) away from the region where 
they are produced and eventually escape into the surroundings 
resulting into an effective decrease in total gas internal energy. 

In presence of cooling, the gas internal energy, (pe)gas, will 
be different from that defined by Eq. (0. Indeed, only con¬ 
tributions due to translational and internal degrees of freedom 
(from H, He and H 2 ) should be included. Conversely, terms cor¬ 
responding to the emission of photons (e.g. ionization, disso¬ 
ciation, roto-vibrational cooling of H 2 molecule) are correctly 
accounted for by the right hand side of Eq. (29 1 in the A term. 
Therefore Eq. ([27]) now becomes 


iprfigas 


(eni + Shi + CHe) 


pksT 

rriH 


(30) 


drii 

dt 


'Y^'KiunjUk-iii 

j-k j 


where, expressions for each of the internal energy components 
(28) ^re given in table [T] A similar contribution to gas internal energy 
(see. Eq. (|29|l) from internal degrees of freedom in presence of 
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radiative losses has also been applied to study the role of molecu¬ 
lar hydrogen in primordial star formation (e.g., Palla et al.|1983 
[Omukai & Nishi|19981 l. 

For the present purpose, we focus only in the chemical evo¬ 
lution of atomic and molecular hydrogen. In particular, the to¬ 
tal hydrogen number density wh includes contributions from 
atomic and molecular hydrogen i.e., hh - rini + 2 nH 2 + nnn- 
Contributions to the electrons density come from ionized 
hydrogen (uhh) and from a small but fixed fraction of metals 
(Z ~ 10 In addition to hydrogen, we also consider helium to 
be present with a fixed mass fraction of 0.027. The mass density 
p - ptitot mp is a conserved quantity whereas the total number 
of particle per unit volume (wtot = + «e) is not as it 

may change due to ionization, recombination and dissociation 
processes. The chemical evolution of molecular, atomic and ion¬ 
ized hydrogen is governed by the equations listed in Table The 
code tracks the formation and destruction of these three species 
based on the temperature-dependent reaction rates specified and 
updates the corresponding number fractions. 


Xhi - 


nm 

Hh 


Xh2 - 


tlH 


Xhii - 


nnii 

nn 


(31) 


In dilute regions such as the solar corona, Eq. ( [28] l can be 
simplified by setting dnildt = 0 as the time scales are such 
that a balance is always maintained between collisional ioniza¬ 
tion and radiative recombination. This condition is known as 
Coronal equilibrium or Collisional Ionization Equilibrium (CIE) 
and it differs from LTE in two aspects: i) it is only valid in di¬ 
lute plasma - unlike the LTE where high density environments 
are required - and ii) the ionization fraction are estimated us¬ 
ing Eq. ( [28] l in steady state and not with Saha fractions. The top 
panel of Eig.|^ shows the concentration fractions as functions of 
temperature obtained by solving the steady state version of Eq. 
( |28l l. The dissociation and ionization of molecular and atomic 
hydrogen, respectively, are clearly evident at the temperatures of 
T ~ 3 X 10^ K and T ~ 10^ K. Such equilibrium values can be 
used to initialize fractions in numerical computations. 

The radiative losses implemented in our model can be writ¬ 
ten as 


2 , for 


A(n, X, T) — Aci -t- Arj^ -F A^otvib "t" ^H2diss ^grain 5 (32) 

where, Aci and Arr are losses due to collisional ionization and 
radiative recombination, respectively (Te§ileanu et al. 2008) 1. 
The remaining terms, Arotvib, AH 2 di.ss and Agrain are associated 
with molecular hydrogen and represent losses due to rotational- 
vibrational cooling, dissociation and gas-grain processes (Smith 
|& Rosen|2003| l. 

The dependence of the cooling function, A(n, X, r)/n^ 
no = 10^ cm~^ is shown in the top and bottom panels of Eig. 

In the top panel, cooling functions are plotted using concentra¬ 
tions values obtained under CIE conditions while the plots in the 
bottom panel corresponds to fixed concentrations obtained for 
T = 4500 K. Eor temperatures r < 10"^ K, the total cooling func¬ 
tion with equilibrium values of hydrogen fractions {black dashed 
curve) is dominated by the contribution of roto-vibrational cool¬ 
ing of molecular hydrogen {blue line). Molecular cooling due to 
dissociation {green curve) and gas-grain processes {red curve) 
have little impact on the total cooling function. For temperatures 
T > 10^ K, most of the hydrogen molecules have dissociated 
into atoms (see Fig|^ and the total cooling curve in Fig. is 
governed by atomic processes like radiative recombination {cyan 
line) and collisional ionization {magenta line). 

The total cooling curve in the bottom panel of Fig. with 
fixed values of hydrogen fractions differs substantially from the 
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panel above. It is dominated by gas grain processes {red line) for 
T < 100 K. For temperatures between 100 K and 10^ K, molec¬ 
ular cooling due to rotational and vibrational processes {blue 
line) plays a vital role. Even at large temperatures, T > 10"^ 
K, the total cooling curve has major contribution from molec¬ 
ular dissociation {green line), for the chosen fixed fraction of 
molecules. The contribution of collisional ionization {magenta 
line) also becomes important at these temperatures. However, 
the cooling due to radiative recombination remain negligible for 
all temperature values due to extremely small and fixed fraction 
of electrons, X/^//. 


3. Numerical Implementation 

In a Godunov-type scheme the MHD equations ([T]l-(|^ are dis¬ 
cretized using a flux-conserving form where the basic building 
block is 

= Uf - — ^ [(A,F,)i,- (A,F,)i_] + Si (33) 

* d 


where U = {p,p\, B, E, pX/^) is our vector of conservative vari¬ 
ables, At" is the time step and Si is a source term. Here we em¬ 
ploy an orthogonal system of coordinates with unit vectors 
(here d - 1 , 2 ,3 or simply x, y, z) where A^i is the area element 
in the d direction, F;/ is the flux computed with a Riemann solver, 
'Vi is the cell volume while i = {i, j, k) is the position of the com¬ 


putational zone in the domain. We remind the reader to Mignone 
et al. P007| l, Te§ileanu et al. ( 2008| l and |Mignone et al.| ( |2012| l 


for more details. Here we focus only on those aspects that cru¬ 
cially depend on the choice of the equation of state, namely: i) 
the computation of the flux via a Riemann solver and ii) the con¬ 
version between internal energy and pressure. 
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Table 2. Summary of the chemistry reaction set. T is the temperature in Kelvin, Tev is the temperature in electron-volts, Tj = T/l x 10^ and T 2 
= T/lOO 


No. 

Reaction 

Rate Coefficient (cm-^s ‘) 

Reference 

1. 

H + e- ^ H+ -(- 2e- 

ki 

= 5.85 X 10-” exp(-157,809.1/T)/(1.0 -t TfO 

1 

2. 

H+ -H e- ^ H + hv 

^2 

= 3.5 X 10-'2(r/300.0)-'’ * 

2 

3. 

H 2 -t — > 2H - 1 - 

h 

= 4.4 X 10-‘°r'’ 35exp(-102,000.0/T) 

3 

4. 

Hj -t H ^ 3H 

ki, 

= 1.067 X 10-‘®r2^‘>‘2(gxp(4.463/Tev)^‘((l- + 0.2472Tev)^-^'^)^‘ 

4 

5. 

Hj -t Hj ^ Hz -t 2H 

^5 

= l.Ox 10-*exp(-84,100/T) 

2 

6. 

H -t- H Hz 

h 

= 3.0 X 10^*’ VT(1.0 -h 0.4 VTz +0.15 + b.2T2 + O.STj) 

5 


References. (1)|C^(|T992) [Eq. 26a]; (2)|Woodall et aL|(|2007ll [UMIST Database] (3)|Galli & Palla|(fT998l( [Eq. H17]; (4)|Abel et al.|p997j [Tab. 
3 Eq. 13]; (5) |Hollenbach & McKee| ( |l979^ [Eq. 3.8] 


10 “^® 



•^H2dis.s 

- Aci 

-^grain 

-^rotvib 

^RR 


10 -“ 



Temperature [K] 



Temperature [K] 


Fig. 3 . Top Panel: Variation with temperature of cooling function aris¬ 
ing from different processes obtained using equilibrium values of hy¬ 
drogen fractions (that vary with temperature as shown in Eig.|^ in units 
of ergs cm^ s“' for a total number density no = 10^ cm“^. Bottom panel: 
Different components of radiative cooling functions with same initial 
number density, no but for fixed values of concentrafion fractions (men¬ 
tioned in the figure). In both panels, the sum of all the components is 
drawn with a black dashed line to obtain the value of A in ergs cm^ s“* 
following Eq. l |^ . 


in a self-consistent way, the initial left and right states through 
wave-curves. For a convex EoS, as defined in Sect. |2.2| the so¬ 
lution to the Riemann problem in hydrodynamics consists of a 
left-facing shock or rarefaction wave, a contact discontinuity and 
a right-facing wave (again either a shock or a rarefaction). No 
compound wave can be formed in the solution. In MHD the so¬ 
lution may involve up to 7 modes: a pair of fast magneto-sonic 
waves, a pair of Alfven waves (or rotational modes) and a pair 
of slow waves and a contact (or tangential) discontinuity in the 
middle. Although exact solutions to the Riemann problem are 
possible, the computational overhead is largely reduced by using 
approximate solvers based on different levels of simplification. 

A first class of solvers heavily relies on characteristic 
(Jacobian) decomposition or computation of Riemann invari¬ 
ants which is strictly connected with the underlying form of the 
conservation law. Typical examples are linearized (Roe-type), 
flux-splitting or two-shock approximate solvers, see the book by 
Toro| (2009]l. Solvers belonging to this class require consider¬ 


able changes when new physical ingredients (such as a different 
EoS) are introduced. In the case of real gases, for instance, ex¬ 
tensions have been presented by|ColeIIa & Glaz ( I985|l;[Glaister 


( |I988b|a| l; [Buffard et al.| ( |200()| l (in the case of the Euler or 
Navier-Stokes equations) while generalization to the MHD case 
have been presented in Dedner & Wesenberg|(|2001 1 . In Fedkiw 


et al. ( I997|l a Roe-type Riemann solver for the solution of ther¬ 


mally ideal, chemically reacting gases has been presented in the 
context of the multi-species Navier-Stokes equations. 

A second class of solvers employs only minimal informa¬ 
tion (typically approximate expressions for the wave speeds) and 
it is based on an application of the integral form of the con¬ 
servation laws which gives a closed-form approximate expres¬ 
sion for the fluxes. Typical examples are the Lax-Friedrichs- 
Rusanov ( Toth & Odstrcil|[l996| l, Harten-Lax-van Leer (HLL, 


see 


Harten et al. 1983 1 solver and their extensions such as HLLC 


dToro et all 1994||Gurski|2004l|Li|2005| ) and HELD ( |Miyoshi &| 


Kusano||20C)51l. Solvers belonging to this class can accomodate 


new changes with minimal efforts by simply changing the def¬ 
inition of the sound speed and, for this reason, will be our pre¬ 
ferred method of choice. Incidentally, we also note that, since 
only an approximate expression for the eigenvalue is needed and 
Fi < 5/3, employing Eq. (|^ with Fi = 5/3 provides an up¬ 
per bound to the actual expression with a only a slight loss of 
accuracy. 


3.1. On the solution of the Riemann Problem 

The Riemann problem is an initial value problem describing the 
decay of a discontinuity separating two constant states. As time 
evolves, the discontinuity breaks into a set of non-interacting el¬ 
ementary waves whose properties are determined by connecting. 


3.2. Conversion between internal energy and pressure 

The computation of the right hand side of Eq. ( [3^ is normally 
carried out using primitive variables customarily defined as V = 
ip, V, B, p). The conversion between U and V requires obtaining 
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Table 3. Equations being solved when converting from pressure to internal energy (p —» pe) or viceversa {pe —> p). Equations marked with 
a * are nonlinear and may require a root-finder approach. Under non-equilibrium conditions the chemical fractions X are independent variables 
available immediately after the hydro update. Under LTE or CIE, X = X(r,p) and this introduces additional nonlinearities in some equations. 


Conversion 

Physical Conditions 

Temperature 

Einal Equation 

(p,X,p) pe 

[Non-equilibrium] 

T = -Kp(X) 

P 

pe = pe(T, X) 

(p,X,pe) p 

[Non-equilibrium] 

pe(T, X) - pe = 0* 


^ Kp(X) 

ip, P) -»pe 

[LTE or CIE] 

T --Kp(T,p) = Q'‘ 

P 

pe=pe(T,p) 

ip,pe) -» p 

[LTE or CIE] 

pe(T,p) - pe = 0* 

P= 

Kp(T,p) 


pressure from total energy density and viceversa. While internal 
energy density is readily obtained from Eq. (|^, the conversion 
p —» pe and its inverse pe —> p strictly depends on the choice of 
the caloric equation of state. 

For the constant-E EoS, these transformations take a small 
fraction of the computational time as the relation between inter¬ 
nal energy and pressure is straightforward and given by 

pe - ^ . (34) 

Note also that the temperature does not explicitly appears in the 
previous definition. 

The situation is different, however, for a more general EoS 
where a closed-form expression between pressure and internal 
energy may not be easy to obtain. From the considerations given 
in the previous sections, in fact, we can write the thermal and 
caloric equations of state as 

pksT 

, ^ mufi{X) ( 35 ) 

e = e(T,X) 

where X may be an independent variable or, in equilibrium con¬ 
ditions, a function of temperature and density. The explicit de¬ 
pendence on the temperature introduces two additional interme¬ 
diate steps, namely; 

1. During the conversion p —» pe one first needs to compute T 
from the thermal EoS (first of Eq. [35] i: 

p m„p(X) 

r =- - -. (36) 

P kg 

Under non-equilibrium conditions, p = p(X) is a known 
function of the gas composition and Eq. ( [36] ) can be solved 
directly. Under LTE or CIE, on the other hand, X -X(T, p) is 
a function of density and temperature and Eq. becomes 
a nonlinear equation in the temperature variable. 

2. During the inverse transformation (pe p), one must first 
solve for the temperature by inverting 

e = e(7’,X), (37) 


3.2.1. Conversion using root finders 

In this approach we employ the exact analytical expressions 
to compute pressure and internal energy as functions 
of temperature or viceversa. For the caloric EoS, numerical in¬ 
version using a root-finder algorithm is required to obtain T from 
pe(T,X) or pe(T,p) as they both are nonlinear functions of the 
temperature. Additionally, under LTE or CIE, the thermal EoS 
must also be inverted numerically to obtain temperature since 
the mean molecular weight introduces nonlinearity. These cases 
are marked with a * in Table |3] 

The root-finder approach results in increased computational 
cost inasmuch the internal energy is an expensive function to 
evaluate. Among different root-solvers not requiring the knowl¬ 
edge of the derivative, we have found Brent’s or Bidders’ meth¬ 
ods to be practical and efficient root-finding algorithms. 

3.2.2. Conversion using Tables 

A second and more efficient strategy can be used when the inter¬ 
nal energy is a function of temperature and density alone (which 
is typically the case under equilibrium conditions, CIE or LTE) 
and it consists of employing pre-computed tables of pressure and 
internal energy, e.g., 

{p)ij = PiTi,Pj), {pe}ij = pe{Ti,pj) (38) 

where i - 0, ■■.,Nc - 1 and j - 0,..., W - 1 are the table indices. 
For convenience, the tables are constructed using equally-spaced 
node values in logT and logp so that log T,/To = iAlogT and 
logpy/po = yAlogp where po and Tq are the lowest density 
and temperature values in the table. Following this approach, 
Eqns. ([36]l-([37| and their inverses are replaced by direct or re¬ 
verse lookup table operations. We point out that, in order to be 
invertible, tables must be obtained from monotone functions of 
their arguments and this is always verified for the EoS of interest. 

When T and p are known, for instance, pressure and inter¬ 
nal energy are obtained by using lookup table followed by two- 
dimensional interpolation between adjacent node values. Thus 
we first locate the table indices i and j by a simple division; 


where, under equilibrium assumptions, the chemical com¬ 
position is a function of temperature and density, i.e., X = 
X(T,p). 


i = floor 


logT/To ) 
AlogT j ’ 


j - floor 


logp/Po \ 
A logp / 


(39) 


Table [^summarizes the relevant equations to be solved. 

In order to cope with the numerical inversion of Eqns. 
and ( |37] i we have considered and implemented two different so¬ 
lution strategies that we describe in the following sections. 


where T and p are the input values at which interpolation is de¬ 
sired. Internal energy (and similarly pressure) is then computed 
as 


pe(T,p) = Sij(x)(l -y) + Sij+i(x)y (40) 
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where x ^ (T - 7’;)/(7’,+i - T,) and i/ = (p - Pj)l(pj+\ - Pj) 
are normalized coordinates between adjacent nodes while S is 


an interpolating spline, see Sect. 3.2.3 

Conversely, when pe and p are known, Eq. ( |40) i must be in¬ 
verted to obtain T. To this end, we first locate the index j using 
the second of ( [39] l with p given. In order to find the column in¬ 
dex i, we first construct the intermediate one-dimensional array 
qi - {pe),j(l-{/)-i-{pe)/_/+ij/with y fixed from the previous search 
whereas i - 0, .■■,Nc- A binary search is then performed on qt to 
obtain i and Eq. ( [40| can be inverted to obtain x. Temperature is 
finally obtained as T = T, -i- x{Ti+\ - T,). The inversion depends 
on the choice of the interpolating function S and it is discussed 
in more detail in the next section. 

The tabulated approach has found to be faster than the gen¬ 
eral root finder method giving considerable speedups, up to a 
factor of 4 for certain problems. 


or Halley’s method which gives cubic convergence; 




2/z(xW)/z'(.rW) 
2[/!'(xW)]^ - /z(xW)/!"(xW) 


(46) 


We have found that both methods hardly never require more than 
3 iterations to achieve an absolute accuracy of 10“". 


4. Numerical benchmarks 

The computational framework presented in this work has been 
implemented in the PLUTO code ( |Mignone et al.|2007[|2012 i. 
Selected numerical benchmarks are now presented with the aim 
of investigating the influence of the EoS on problems of astro- 
physical relevance as well as comparing the computational effi¬ 
ciency of the proposed numerical approaches. 


3.2.3. On the choice of the interpolant 

Eor linear interpolation, we simply use 


Si,j(,x) = {pe),+i,yX; -H {pe),j(l - X/) 


(41) 


in Eq. (40i so that internal energy and pressure become bilinear 


functions of temperature and density. When pe is given, on the 
other hand, the inverse function is required to compute tempera¬ 
ture and one needs to solve Eq. (40 1 with >S(x) given by Eq. (41 1 
to obtain x. Since the resulting equation is linear in x one readily 
finds: 


/ - fi.j - - fi.j) 


l/(/i-t-lj+l fij+^ fij^ fi+^J fij 


(42) 


where / = pe, fij = {pe)/,y or / = p, fij = {p}ij. Despite its 
simplicity, however, bilinear interpolation may generate thermo¬ 
dynamically inconsistent results since the positivity of the fun¬ 
damental derivative (Eq. [T6] l can be violated at node points in 
the table where derivatives are not continuous. This generates 
composite waves in the solution due to a loss of convexity of the 
caloric EoS as discussed in Sect. |2.2| Indeed, we have verified 
that this pathology worsen as the number of points in the temper¬ 
ature grid is reduced. An illustration of the problem is reported 
in Appendix [B] 

To overcome this limitation, we have also implemented a cu¬ 
bic spline when interpolating along the temperature grid so that 

<Sij(x) - Oijx^ + bijx^ + Cijx + dij , (43) 


where the coefficients a,b,c and d are computed by ensuring 
that the cubic is strictly monotonic in the interval [T,, r,+i] see 
Appendix The inverse function requires finding the root of 
a cubic equation on a specific interval and, since the spline is 
monotone by construction, only one root is always guaranteed 
to exist. Thus combining Eq. ( |43| ) with Eq. ( |40] l one obtains the 
following equation; 

h(x) - dijx^ + bijx^ + Cijx + dij - f - Q (44) 


where dij - a,,/! - ^) -i- aij^iy and similarly for the remaining 
coefficients. Although standard analytical solvers (e.g. Cardano) 
may be used, we have found root-finder methods to be more 
robust, accurate and computational efficient for the purpose. 
Starting from an initial guess given by Eq. ( |42| i, one can either 
use Newton-Raphson method to achieve quadratic convergence. 


[k+\] ^ m _ hix^’^'^) 
/!'(xW) ’ 


(45) 


4.1. Sod shock tube 

In our first test, we consider the standard Sod shock tube on the 
unit interval x e [0,1] with a uniform grid resolution Ax = 10“^. 
The initial setup consists of two fluids initially separated by a 
discontinuity located at x = 0.5. Density and pressure in the two 
regions are given hy (p, p) - (1,1) for x < 0.5 and (p, p) - 
(1/8, 1/10) for X > 0.5 while the velocity is zero everywhere. 
We choose our physical units in such a way that the reference 
velocity is vq - 5.25 x 10^ Km/s so that the temperature, in 
Kelvin, is expressed by T = {plp)v^pmulkB. With this choice, 
the temperatures in the left and right regions become Ti ~ 3852 
K and Tr ~ 3345 K. The final solution at f = 0.2 is shown in 
Fig. 12 for the adiabatic run (i.e. without explicit cooling). The 
test is repeated also by including thermal losses and the solution 
is plotted in Eig. |2 

4.1.1. Acdiabatic Shock Tube 

Eig.|2compares the spatial variation of density, pressure and ve¬ 
locity for an ideal EoS {solid lines) with that of a H/He EoS 
dashed lines) without explicit cooling. The solution comprises, 
from left ro right, a rarefaction wave, contact discontinuity and a 
right-going shock. Pressure and velocity are always continuous 
across the contact wave. In the case of the gas with H/He EoS we 
obtain a larger compression ratio (a; 2.8 compared to a; 1.84 of 
the ideal gas) and the shock propagates slightly slower. Similarly 
the tail and the head of the rarefaction wave both propagate at 
a smaller speed and this owes to a reduced value of the sound 
speed. The density jump across the contact wave is largely re¬ 
duced with the H/He EoS. 

It is important to understand that in the case of an ideal EoS, 
only translation degrees of freedom contribute to the internal en¬ 
ergy which is a linear function of the temperature, see Eq. Gl- 
In the case of the H/He EoS, however, the internal energy also 
takes contributions from additional degrees of freedom due to 
rotations and vibrations of di-atomic molecules (like H 2 at low 
temperatures). These will not correspond to an increase in tem¬ 
perature. Therefore, across a shock wave, the upstream kinetic 
energy will become available to the system not only to raise 
the temperature but also to populate the vibrational and rota¬ 
tional levels of the molecules (at lower temperatures), dissociate 
molecules and eventually ionize atoms (at larger temperatures). 
Hence we expect the overall increase in temperature to be re¬ 
duced in the case of H/He EoS as compared to that of an ideal 
EoS. 
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Fig. 4. The figure shows the variation of density p (red). Pressure P 
(green), and velocity (blue) along the X-axis (in code units) for a 
standard Sod Tube test (without explicit cooling) at time t = 0.2. The 
values obtained using an ideal EoS are shown as solid lines while that 
obtained using a GammaLaw EoS are shown as dashed lines. 


Finally we point out that while the employment of an ideal 
EoS gives a scale-free and self-similar solution to the Sod 
shock tube, the same argument does not hold for the H/He EoS 
which implicitly contains temperature scales corresponding to 
the above mentioned transitions. 


4.1.2. Radiative Shock Tube 


We have repeated the Sod shock tube including the explicit non¬ 
equilibrium cooling described in Sect. 2.6 We consider two ini¬ 
tial conditions corresponding to different values of temperature. 
In the first case the temperature of the left and right states is 
set X.O Ti - 400 K and Tr - 200 K, respectively. The second 
case corresponds to hotter gas, Ti — 4500 K and Tr - 3500 K. 
Molecular and atomic hydrogen fraction are initially set to their 
equilibrium values at the local temperature (see Eig.|^. 

Fig.|5] shows the density and temperature (T/p, p being the 
mean molecular weight). In all panels, solution obtained using 
an ideal EoS is shown as a black line, whereas that obtained 
using H/He EoS with a red line. 

Eor smaller initial temperatures (top panels), the density and 
temperature jumps differ by a modest factor (=» 20%) and the 
positions of the waves is essentially the same. The differences 
are not as large as in the adiabatic case since the H/He EoS, now 
given by Eq. ( |30| l, only includes terms related to the internal de¬ 
grees of freedom, see Sect. 2.6 Terms related to dissociation and 
ionization correspond to energy losses and are thus accounted 
for by the cooling function A. These terms are common to both 
EoS. The presence of molecules at low temperature gives a con¬ 
siderable contribution to the internal energy budget. Hence part 
of the thermal energy available at the shock serves in populating 
the internal molecular degrees of freedom rather than raising the 
temperature. As a consequence, radiative losses are slightly less 
efficient for the H/He EoS. 




Fig. 5. The figure shows variation of density (left panels) and temper¬ 
ature (right panels) at the final stage of Sod shock tube test for cases 
which have included explicit cooling. The black line represent values 
obtained from ideal EoS, while those obtained using H/He EoS are 
shown in red. The top and bottom panels differ in the value of their 
initial temperature on either sides of the interface at x = 0.5. 


For larger initial temperatures (bottom panels) differences 
are negligible and this behavior owes to a reduced cooling length 
scale. Indeed, in presence of a shock wave, the gas behind the 
front undergoes a rapid increase of the temperature while the 
level populations of H 2 do not occur instantaneously ( |Flower| 
et al.|20d3 1. Since the cooling time scale is shorter than the av¬ 


erage time scale over which level population takes place, inter¬ 
nal energy changes are dominated by radiative losses. Since the 
cooling function A is the same for both EoS the final solutions 
are very similar. 

These different test cases elucidates the importance of treat¬ 
ing temperature-dependent EoS both in case when the system 
is in equilibrium and also when explicit radiative cooling is in¬ 
volved. 


4.2. MHD Shock Tube 


Next we consider the collision between two magnetized fluids 
moving in opposite directions. The problem is studied on the 
unit interval x e [0,1] with initial condition given as follows: 


(p, Vjc, Bx, 


(5, 5, 5, 8, 10) for 
(1, -15, 5, 3, 1) for 


.r < 0.5 
.r > 0.5 


(47) 


The fluid on the left has larger inertia than the fluid on the 
right and it is also more magnetized. We solve the problem by 
using linear interpolation on characteristic variables, the HELD 
Riemann solver and an adaptive resolution consisting of a base 
grid of 128 zones and 7 levels of refinement with consecutive 
jump ratios of 2 except for refining level 3 where the grid jump 
is 4. 

The solution is shown in Fig. at f = 0.08 for the ideal 
and H/He EoS. In order to carry out a fair comparison, we have 
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Fig. 6. The figure has 8 panels, going from top left to bottom right the panels show solution at time t = 0.08 of density, pressure, internal 
energy density (p e), temperature, X and Y components of velocity, Y component of magnetic field and adiabatic index respectively. In each panel, 
quantities obtained for ideal EoS are shown in (black) while the red line denotes values obtained using H/He EoS. 


chosen the adiabatic index of the ideal EoS to be E = 1.447 so 
that the gas internal energy has the same value for both EoS at 
f = 0. The wave pattern consists of a pair of outer fast magneto- 
sonic shocks enclosing two slow magneto-sonic shocks and a 
contact wave in the middle. 

Overall, shocks are stronger and propagate faster for the 
ideal gas case leading to the formation of a more extended 
Riemann fan owing to a larger value of the sound speed (Eq. 
[T0| i. This is evident by looking at the profiles of the first adia¬ 
batic index shown in the bottom right panel of Eig. For the 
H/He EoS, Ti drops to « 1.04 in the shocked regions while it 
remains constant for the ideal EoS. This effect becomes more 
pronounced at the two slow shocks which, in the case of H/He 
gas, propagate to the left at very small velocities forming a very 
thin shell of compressed material with large internal energy. 

In spite of the reduced shock strength, however, the density 
compression attained at shocks is nearly the same or even larger 
for the H/He gas. This behavior can be understood by notic¬ 
ing that lower values of the specific heat ratio should support 
larger density jumps. Indeed, the effective T-value, defined as 
Teff = plipe) + 1 (dashed line in the bottom right panel of Fig. 
1^, becomes smaller with the H/He EoS thus allowing for com¬ 
parable or larger compression ratios (in density and magnetic 
field) even if the shock is weaker. 

The EoS largely determines how the upstream kinetic energy 
is being converted into internal energy and heat: for the ideal 
gas an increase in internal energy corresponds to an increment 
in pressure and temperature while, in the range of temperature 
considered here (10^ K < T < 10"^ K), the same does not hold 
for the gas with H/He EoS. In this case, in fact, the downstream 
internal energy is employed to excite internal vibrational and ro¬ 
tational levels without a corresponding increase in pressure or 
temperature. 


4.3. Blast Wave 

We now consider a two-dimensional blast wave in Cartesian co¬ 
ordinates on the square domain -Lo/2 < x,y < Lo/2 where 
Lo = 2.5 X 10'^ cm. An over-dense and over-pressurized region 
is initialized at the center of the domain (x - y - 0) inside a 
circular region of radius ro/Lo = 0.1 where density and pres¬ 
sure are set respectively to pipa - 10 and pjpa - 20 in units 
of the corresponding ambient values, pa - tiamn and pa - paV^ 
(here we use tia - 10^ cm“^ and uq = 2.5 km/s). The compu¬ 
tation is carried using adaptive mesh refinement (AMR) using a 
base grid of 256 x 256 grid zones and 4 levels of refinement (ef¬ 
fective resolution of 4096^ zones). The MUSCL-Hancock with 
piecewise linear interpolation and the HLLC Riemann solver are 
employed. 

The solution features an outermost forward shock wave fol¬ 
lowed by a contact discontinuity and a reverse shock as shown in 
the top panel of Fig.[7]at t = 0.15. In analogy with the Sod shock 
tube problem, we observe that the size of the spherical blast wave 
is smaller when the H/He EoS is employed although the shock 
strengths (measured in terms of pressure jumps) are comparable 
for both EoS. The compression ratio across the forward shock 
is noticeably larger for the H/He EoS (~ 2.7) than for the ideal 
EoS (~ 1.7) owing to a reduced value of the equivalent adiabatic 
index. 

An important difference lies in the structure of the contact 
wave as it can be noticed in the closeups in the bottom panel of 
Fig. □ In the ideal EoS case, the density contrast across the con¬ 
tact wave is favorable to the onset of the Rayleigh-Taylor insta¬ 
bility while the same structure is stable for the H/He EoS since 
the density contrast is reversed (heavy fluid in front, light fluid 
behind). Thus a shell of high density material is swept between 
the outermost shock and the contact wave. 

From the computational perspective, we have compared the 
ideal and H/He EoS in terms of speed. Our results show that the 
average wall clock time per numerical step is 0.05 s with the 
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Fig. 7. Top panels: Comparison of density p (in code units) for a hydrodynamic, spherical blast wave at time f = 0.15. The left panel shows the 
final shock structure with an ideal EoS. The density obtained using the H/He EoS is shown in the top right panel. Alongside each panel a zoomed 
in view of density around the contact discontinuity is shown. The ideal EoS shows presence of Rayleigh-Taylor instability while its completely 
absent in case of H/He EoS. The bottom panels from left to right compare ID cuts at the mid-plane of density p, pressure p, temperature (T/p) and 
velocity (Vx) for ideal EoS (black line) and H/He EoS (red line). 


ideal EoS while we found 0.25 s and 0.9 s for the tabulated and 
the root finder approaches, respectively, in the case of H/He EoS. 
As expected the computation with constant E are the fastest ones. 
However, it is interesting to note that pre-computed tabulated 
values give a speed up of about 4 times than that obtained using 
the root finder approach, while maintaining essentially the same 
accuracy in the final solution. 


4.4. One-dimensional pulsed molecular jets 

As an astrophysical applications, we consider in the next test the 
propagation of velocity pulsations in a ID molecular jet model 
that includes radiative losses. Multidimensional extensions of 
this study will be considered in forthcoming papers. 

In order to resolve the thin post-shocked regions we em¬ 
ploy adaptive mesh refinement to enhance resolution in high- 
temperature and dense regions. The computational domain ex¬ 
tends from z = 0 up to z = 32 in units of the jet radius rj ~ 75 
AU and it is covered by a base grid of 256 grid zones. We em¬ 
ploy 10 levels of refinement yielding an equivalent resolution of 
262,144 zones. The domain is initially filled with a fully molec¬ 
ular (Xh 2 ~ 0.5) ambient representing a young star-forming core 
at the temperature of Ta = 50 K and decreasing density so that 


Paiz) = 


Pa 

1 + fe/zo)^ ’ 


Paiz) = 


pksTg 

prriH 


(48) 


In the previous equation pa is the ambient density at z = 0 while 
Zo = lOrj. An over-dense jet (pjlpa - 3) is injected from the 
nozzle at z = 0 with temperature Tj = 10^ K resulting in an 
over-pressurized jet (pj !^ 60). The initial hydrogen fractions 
at the nozzle are estimated based on equilibrium shown in Eig.lS 
while the jet density is assumed to be pj - rijinn with rij - 1(T 
cm“^. The injection velocity has sinusoidal perturbations of the 


form Vj = Co (1 + 0.25 sm(2ntlTp)), with base velocity uq = 70 
km s ' and pulsation period Tp = 40 years. 

As the system evolves, velocity pulsations steepen into a 
chain of forward/reverse shock pairs, the first of which becomes 
the strongest one reaching temperatures of T ~ 10^ K where hy¬ 
drogen is quickly ionized. Under these physical conditions the 
two EoS yields essentially the same results since thermal kinetic 
motion is characterized only by translational degrees of freedom. 

Successive pulses, however, form weaker shocks with tem¬ 
peratures of few thousand K, see Eig. As the cooling time 
is much smaller than the dynamical time, the material behind 
the shocks start to cool radiatively in a very efficient way. 
Temperature peaks immediately downstream of the forward and 
reverse shocks and quickly drops afterwards forming extremely 
thin layers of hot material (see the enlargement in the right panel 
of Eig.|^. The layer of fluid between the forward/reverse shock 
pairs becomes nearly isobaric with a high density pile. We point 
out that the numerical resolution of such layers is crucial in or¬ 
der to correctly capture the physical processes (such as ioniza¬ 
tion and radiative losses) taking place on different time scales. 
The importance of grid resolution has also been addressed by 
Tegileanu et al.| ( |2009| l in the context of time-dependent radia¬ 
tive shocks subject to atomic cooling. 

As it has been elucidated in earlier tests, the transforma¬ 
tion of upstream kinetic energy into internal energy depends on 
the EoS considered. While pressure is essentially the same for 
both EoS, the temperature values reached immediately behind 
the shocks are larger for the ideal EoS gas than for the H-He 
EoS because part of the internal energy is employed to dissoci¬ 
ate molecules. As a consequence the amount of atomic hydrogen 
produced during the dissociation process is halved (~ 2.7% for 
the H-He EoS vs ~ 5.9% with the ideal EoS). 
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Fig. 8. Results from the test studying the propagation of supersonic pulsations in molecular medium. Top panel shows the comparison of density 
(in code units) for ideal EoS {black line) and H/He EoS {red line) on a log scale at time ? = 192 yrs. The bottom panels compares zoomed up 
quantities for a particular knot marked within the rectangle in the top panel using the same color keys. These quantities from left to right are 
density p, pressure p, HI fraction X^i and temperature T/p. A zoomed inset of the forward shock obtained from H/He EoS is shown as well with 
red dots to demonstrate the importance of resolving the shock (see text). 


4.5. Axisymmetric MHD jets 


In this test problem we consider the effect of the H/He EoS on 
the propagation of an axisymmetric MHD jet using 2D cylin¬ 
drical coordinates. For simplicity, we consider the MHD equa¬ 
tions without the radiative loss term. The computational domain 
is defined by r 6 [0,10] and z e [0,60] and discretized using 
256 X 1536 grid zones. It is initially hlled with a static ambient 
medium with constant density pa - ninria (na = 10^ cm“^) and 
pressure pa = (3/5)paUy where i>o = 1 km/s is our unit velocity. 
The ambient medium is threaded by a constant vertical magnetic 
held Bj. = ^J2cr^pa where cTj. = 1/4 is the magnetization param¬ 
eter corresponding to the vertical held component. 

A supersonic inhow is imposed through a circular nozzle at 
the lower boundary (z - 0) with radius Rj = 20 AU where the 
beam is inject with the same density as the ambient medium and 
is injected with a constant vertical velocity = 150 km s“'. 
The jet carries an azimuthal magnetic held with strength given 
by ( |Te§ileanu et al.|20()8| 


B^R) 


R 

Bm— for R < a 
a 

a 

Bm — otherwise, 
R 


(49) 


where = ^Acr^Pal{a^(t -41na)), cr^ - I controls the 
strength of the azimuthal magnetic held and a - 0.8 is the mag¬ 
netization radius. Radial balance across the jet beam leads to the 
following (thermal) pressure prohle: 


p{R) ^ Pa +bI, 


1 - min 



(50) 


where the integration constant is chosen in such a way to have 
pressure balance across the jet border. 


As the jet enters into the ambient medium, it immediately 
forms a bow-shock that pushes the ambient material to its sides. 
The overall morphology of the jet is determined primarily by the 
Mach number, the density contrast (a result that has also been 
conhrmed in laboratory experiments ( |Belan et al.||2013| |2014| )) 
and the magnetic held strength. The processed material gets 
heated and forms the cocoon as shown in Fig.|^at t ~ 65 yrs for 
the ideal (left) and H/He EoS (right). The left half of each panel 
shows the logarithmic values of temperature in Kelvin, while the 
right panel shows the corresponding value of density in units of 
the ambient density. For both EoS considered here, the largest 
temperature is achieved at the bow-shock where ~ 2.9x 10^ 
K for the ideal EoS while 2 x 10^ K for the H/He EoS. 

On the contrary, the maximum density obtained with H/He EoS 
is Hmax ~ 1.8 X 10^ cm“^, larger than its ideal counterpart by a 
factor of ai 10. 

The two simulations differ mostly in the overall morphol¬ 
ogy. While the ideal EoS gas tends to produce a larger cocoon, 
the H/He gas results in a thinner and colder cocoon, a struc¬ 
ture that resembles that of a radiative jet. Once again, we notice 
that the shocked ambient medium contributes more to raising the 
internal energy rather than the temperature inside the cocoon. 
This results in a lowered thermal pressure to support the mate¬ 
rial feeding the cocoon leading to the formation of a thin dense 
layer on its sides. This thin layer also shows formation of rolled 
up Kelvin-Helmholtz vortices. 

Similar to the blast wave test (see |4.3| l, we have performed a 
computational efficiency comparison. We found the average wall 
clock time per step to be 0.048 s for the constant-E EoS while, 
for the H/He EoS, the timing were 0.095 s using tabulated ap¬ 
proach and three times as much using the root-hnder conversion 
method (Afav = 0.29 s). These results confirm the trend already 
established for previous tests. 
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Fig. 9. Comparison of logarithmic values both of density, p (right half in blues) and temperature in Kelvin (left half in red) for a 2D axisymmetric 
jet without explicit cooling. The left panel shows the jet structure with an ideal EoS, while, that obtained from a H/He EoS is shown in the right 
panel. 


5. Summary 

In this work, we have revisited fundamental thermodynamic 
properties for the modeling a gas mixture in terms of its thermal 
and caloric equations of state. This has been achieved by consis¬ 
tently including temperature-dependent physical processes such 
as ionization, dissociation and level population of internal de¬ 
grees of freedom. This approach has shown to considerably im¬ 
prove over the (widespreadly used) constant-T EoS in those as- 
trophysical environments where additional physical processes 
play a vital role. In this respect, we consider the H-He EoS that 
takes into account various atomic and molecular processes in 
the context of equilibrium conditions as well as in the presence 
of non-equilibrium cooling. 

The paper also presents a detailed numerical implementation 
of thermally ideal and general caloric convex EoS for Godunov- 
type numerical schemes. In particular, it has been shown that 
conservative schemes require two major modihcations. The hrst 


one concerns the Riemann solver and only Jacobian-free solvers 
have been considered in the present context. Eor such schemes 
only a change to the sound speed is necessary thus making their 
extension to a more complex (convex) EoS straightforward. The 
second modihcation relates to the conversion between total con¬ 
served energy and gas pressure. Since a general thermal EoS can 
now be nonlinear functions of the temperature, a closed-form 
expression between internal energy and pressure cannot be ob¬ 
tained and the problem must be treated numerically. To cope 
with this, we have explored two alternative strategies: the hrst 
one relies on numerical inversion through root-hnder methods 
whereas the second one is based on a combination of lookup 
table and interpolation. We found the hrst approach to be accu¬ 
rate but also computationally more expensive (by a factor ~ 20) 
when compared to the Ideal EoS. The second approach, best 
suited for equilibrium conditions or tabulated EoS, yields sim¬ 
ilar accuracy and it largely reduces the computational workload 
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(a factor 3 ~ 4 in the presented benchmarks). Cubic spline inter¬ 
polation (for the tabulated approach) is preferred over bilinear 
interpolation in order to ensure thermodynamic consistency so 
as to avoid the loss of convexity and the ensuing formation of 
composite waves in the solution of the Riemann problem. 

The proposed numerical benchmarks indicate that the differ¬ 
ences between the ideal and H-He gas can be considerable when 
the gas temperature lies in a range favorable to dissociation or 
ionization processes as well as to populate levels corresponding 
to roto-vibrational degrees of freedom. The differences are par¬ 
ticularly evident in problems involving shock propagation. For a 
mono-atomic ideal EoS (constant F), in conditions of local ther¬ 
modynamic equilibrium, the energy generated during the shock 
impact can only be distributed into translational degrees of free¬ 
dom thus corresponding to a temperature increment. Conversely, 
for the H-He EoS, since additional degrees of freedom are avail¬ 
able, the internal energy of the post-shock gas becomes available 
for populating rotational and vibrational levels rather than result¬ 
ing in an increase in temperature. In addition, the employment of 
such an EoS for temperatures 10^ < T <10^ (depending on the 
density) results in an lower effective adiabatic index and larger 
compression ratios are observed at shocks. 

These features significantly alter the structure of the solution 
of the Riemann problem in terms of wave strengths and posi¬ 
tions and frequently result in the formation of thin-shell of high 
density gas. In a 2D environment, for example, the density con¬ 
trast across the contact wave may even reverse when switching 
from the ideal to the H-He EoS thereby completely quenching 
the growth of the Rayleigh-Taylor instability (see the blast wave 
test). 

Under non-equilibrium conditions when optically thin cool¬ 
ing is introduced, a similar trend is observed although differ¬ 
ences between the EoS are visible mainly at lower temperature 
range {T ~ 100 K). In such a case, terms related to the molecu¬ 
lar degrees of freedom are added to the internal energy while the 
terms responsible for radiative cooling are accounted for by the 
cooling function. We stress out that very high spatial resolution 
may be required to resolve the post-shock regions since the cool¬ 
ing length becomes much smaller than the overall spatial scale. 
On the contrary, in the high temperature regime where hydrogen 
is mostly ionized, the behaviors of the ideal and H-He EoS be¬ 
come similar as the expressions for the internal energy (and the 
radiative losses) are essentially for the same. 

The numerical methods described here have been imple¬ 
mented in the PLUTO code ( [Mignone et al.|20W i starting with 
version 4.2 and are now publicly available. We believe that this 
set of tools can become useful in the investigation of various as- 
trophysical systems like planet formation, molecular jets from 
young stars, gravitational collapse of molecular cloud, galaxy 
formation due to cold flows etc. 
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Appendix A: Fundamental Gas Derivative 


Using standard thermodynamic relations, equation (16i can be 


expressed as a function of p and temperature T (Nannan et al. 
|20l3l ) 


@ - 


2clp- 


-(Fi H- F2 -H F3) 


(A.l) 


where the thermodynamic speed of sound, Cj, is computed from 


dP\ T IdPV 
dpjj C„p2 (57 


(A.2) 
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and the three F’s are given by, 


ri 



+ 2p^ 



ra 




'c«W)p 


d^P \ 



1 IdP 


T \dT 


I dT IJ 


(A.3) 

where, C„ denotes specihc heat capacity at constant volume. 


Appendix B: Spurious formation of composite 
waves 


It has been shown in Section 3.2.2 that lookup table methods 
can be faster than nonlinear root-hnder methods (Sect. |3.2.1| i 
although care must be taken to ensure consistency with ther¬ 
modynamic principles. Indeed, as stated in Sect. 3.2.3 inaccu¬ 


rate interpolation may give rise to a local loss of convexity with 
the ensuing formation of composite waves in the solution of the 
Riemann problem. This failure can be ascribed to the choice of 
a low-order interpolant and can be quantihed by numerical com¬ 
putation of the fundamental gas derivative, Eq. 

Here we present the analysis for a partially ionized hydrogen 
gas in LTE equilibrium with internal energy given by 


e = -kbT +x 


(B.l) 


where AT/ = Tj+i - Ti and, for simplicity, we have dropped 
the index j spanning along the density grid. Cubic polynomi¬ 
als are frequently used as they allow to smoothly ht a set of 
data values with specihed derivatives at each endpoint. Although 
continuity of the hrst and second derivatives across all polyno¬ 
mial segments are normally invoked when deriving the coeffi¬ 
cients of Eq. Q, we shall pursue a different appr oach where 


continuity is sacrihced in favor of monotonicity (Wolberg & 


Alfy|2002]l. The resulting spline will be the smoothest curve that 


passes through the control points while preserving monotonicity 
of the dataset but not necessarily continuity of the second deriva¬ 
tive. It can then be shown that the coefficients are 


where 


a,- = ATi {-Inti + fI + 
bi = ATi{3mi-2f' - f'^^) 

ci = atj; 

di = fi 


Mi = 


/i+i - f i 
Tm - Ti ■ 


(C.2) 


(C.3) 


The cubic is monotone if there is no change in sign in the 
interval [xi,Xi+i]. If this holds, then a necessary conditions is 
that 


sgn(/;') = sgn(/;.'^j) = sgn(m,) (C.4) 


Dehning a,- = flI'n, andjS, = fl^^firii, and assuming that m, 4 0, 
the sufficient conditions for monotonicity can be summarized as 
follows: 


where x is the ionization energy for hydrogen. 

The internal energy is characterized by a rapid increase 
around T ^ 4000 K (corresponding to ionization) and it is plot¬ 
ted in the left panel of Eig. |B.1| together with piecewise polyno¬ 
mial approximations using linear (Eq. red curve) and cubic 
(Eq. 43 green curve) splines. The table was constructed by sam¬ 
pling the temperature range [1,10^] K with Nc = 350 points 
using a uniform log spacing A log T - 8/350. Clearly, the tran¬ 
sition region is poorly sampled and the linear interpolation man¬ 
ifestly reveal that the hrst derivative is not continuous. This be¬ 
havior has dramatic consequences when plotting the prohle of 
the fundamental derivative computed using central differences 
and shown in the middle panel of Eig. B. 1 Large spikes are ev¬ 
ident at the juncture points between contiguous piecewise linear 
polynomials while this effect is not seen using monotone spline 
interpolation. The large over- and under-shoots fail to preserve 
local convexity leading to negative values of 0. When applied 
to a shock-tube test problem, this loss of convexity results into 
composite waves as observed in the rarefaction branch shown 
in rightmost panel of the hgure. The number of spikes can be 
reduced by increasing the number of sample points for linear in¬ 
terpolation to 512 although it is completely absent when using 
cubic spline interpolation. 


1. Eq. ( |C.4|l and o',- H- /?, - 2 < 0; 

2. Eq. ( |C.4) and a, -i- 2jS, - 2 > 0 and one of the following 
conditions: 

a) 2ai H- ySi - 3 < 0; 

b) ai + 2j3i - 3 < 0; 

c) af + aiiPi - 6) -H {Pi - 2>f < 0. 


If monotonicity cannot be satished inside the interval or the 
spline does not have enough curvature, we simply revert to 
linear interpolation. A test for the latter condition is simply 
\d^fi(T)ldT^\ < e\fi(T)\ computed at the interval midpoint T = 
(Ti + Ti+i)f2 and e = 10“*® is a small number. 


Appendix C: Monotone Spline 

Given a set of values fi - fiiTi) and theirs derivatives fil — 
dfiIdx{Ti) dehned at node points Ti we construct a piecewise 
cubic approximation on each interval [T,, r,+i] in the form 


fii{T) 


T-Tj 

ATi 


+ bi 


T -Ti 
ATi 


+ Ci 


T-Ti 

ATi 


+ di 


(C.l) 
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Fig.B.l. Comparisons between exact and tabulated approaches for the partially hydrogen gas EoS. The black dashed line corresponds to compu¬ 
tations obtained with the iterative root-finder approach outlined in Section [3.2.l| while the red and green solid lines corresponds to the tabulated 
approach presented in Section [3.2.2| Left: gas internal energy (Eg. |B.1| ( around the transition region. Middle’, fundamental gas derivative @ (Eg. 
|A.l^ computed using standard central dilferencing. The large spikes m the red curve are due to discontinuous first derivatives using linear interpo¬ 
lation between adjacent node values in the table. Right: Density profile in the Sod shock tube at f = 0.2 using the exact and spline approaches. It 
is clear the appearance of spurious compound structures in the rarefaction wave when a linear spline is used to approximate the table. 
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